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CYCLIC MAGNETIC ACTIVITY DUE TO TURBULENT CONVECTION IN SPHERICAL WEDGE GEOMETRY 
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ABSTRACT 

We report on a simulation of turbulent, rotating, stratified, magnetohydrodynamic convection in spherical 
wedge geometry. An initially small-scale, random, weak amplitude magnetic field is amplified by several 
orders of magnitude in the course of the simulation to form coherent large-scale fields in the saturated state of 
the dynamo. The differential rotation is solar-like (fast equator), but neither coherent meridional circulation 
nor near- surface shear layer develops in this run. Regions of strong magnetic field propagate towards the poles 
at high latitudes and towards the equator at low latitudes, reminiscent of the solar cycle. 
Subject headings: Magnetohydrodynamic s - convection - turbulence 



1. INTRODUCTION 

The solar magnetic field exhibits a quasi-periodic cycle 
with a period of approximately 22 years. This cycle is man- 
ifested by the appearance of sunspots in low latitude activ- 
ity belts that migrate towards the equator as the sunspot cy- 
cle progresses. Reproducing this behaviour remains a major 
challenge to theoreticians. Mean-field mod els, where small- 
scale turbulent effects are parametrized (e.g. lKrause & R adler 
1980), have reproduced many aspects of the solar cycle, but 
with broadly varyi ng assumptions for the variou s parametriza- 
tions (see, e.g. Dikpati & Charbonneau 1999; Ossendriiver 
l20Q3l:lKapvla etal . 2006; K itchatinov & Qlem skov 20ll|). 

Another, computationally much more demanding, but phys- 
ically more consistent route is to solve the equations of mag- 
netohydrodynamics directly without resorting to ill-defined 
parametrizations for the small scales. In practise, how- 
ever, realistic Reynolds and Rayleigh numbers, describing 
the effects of molecular diffusion with r espect to advection , 
are not accessible to simulations (e. g. [Chan & Sofia 1986; 
Miesch & Toomre 2009; Kapyla 201 1). The usual approach is 
to enhance the diffusion coefficients to levels that are compu- 
tationally feasible while striving to maximize the resolution. 

Early spherical shell simulations were able to produce a 
solar-like rotation profile, i.e. one with ' 'eq ua torward acc el- 
eration," and osc illatory large-scale dynamos (Gilman 1983; 
GlatzmaieU 11985b . However, the direction of propagation of 
the dynamo wave was towards the poles, in contradiction to 
the Sun. This can be qualitatively explained by a Parker 
dynamo wave with positive radial shear near the equator in 
conjunction with negative kinetic h elicity, or a positive a- 
effect, in the northern hemisphere (Parker HT955l) . More so- 
phisticated simulations with solar rotation rate and luminos- 
ity have failed to produce strong large-scale mag netic fields 
(Bru n et all 120041) or clear cyclic behaviour dMiesch et al.l 
2011]). These runs omitted a stable layer below the convection 
zone. When such a layer is added, non-oscilla tory large-scale 
fields are found also for solar parameters (Browning et al.l 
2006). Later, oscillatory solutions have been obtained for 
more realistic boundary c onditions and subgrid- scale model- 
ing dGhizaru et al.ll2010t iRacine et~aT1l201 ll) . These are the 
most solar-like solutions so far, but also in their model the 
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activity is at too high latitudes and the activity belts do not 
propagate towards the equator. When the rotation rate is in- 
creased from the solar value in runs without an overshoot 
layer, first stable w reaths of strong large-scale fields appear 
( Bro wn etaLl 12008b , and at even more rapid rotation , pole- 
ward migrating ac tivity is found dBrown et al.l 120 lOi 120111; 
iKapvla et al]l20lQh . 

We report here results from a simulation of turbulent con- 
vection in spherical wedge geometry with solar-like equa- 
torward acceleration that exhibit, for the first time, equator- 
ward migrating magnetic activity near the equator and a polar 
branch at high latitudes. Th e numerical model is the same as 
that in Kapyla ^eTaT] J2011al) but here we add magnetic fields. 

2. THE MODEL 

We model a segment of a star, i.e. a "wedge", in spherical 
polar coordinates, where (r, <p) denote radius, colatitude, 
and longitude. The radial, latitudinal, and longitudinal extents 
of the wedge are 0.7R < r < R, < < ir - 6> , and 
< <j) < 0o, respectively, where R is the radius of the star. 
Here we take Oq = tt/8 and (/>o = tt/2. 

We solve the compressible hydromagnetics equations, 
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T— = - [V • (if VT) + V • (xtpTVs) + 2vS 2 ] , (4) 

where A is the magnetic vector potential, u is the velocity, 
B = V x A is the magnetic field, J = /1q 1 V x B is the cur- 
rent density, r\ is the magnetic diffusivity, /in, is the vacuum 
permeability, D/Dt = d/dt + u • V is the advective time 
derivative, v is the kinematic viscosity, K is the radiative heat 
conductivity, xt is me unresolved turbulent heat conductivity. 
p is the density, s is the specific entropy, T is the tempera- 
ture, and p is the pressure. The fluid obeys the ideal gas law 
with p = (7 — l)pe, where 7 = cp/cy = 5/3 is the ratio of 
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specific heats at constant pressure and volume, respectively, 
and e = cyT is the internal energy. The gravitational ac- 
celeration is g = —GMr/r 2 , where G is the gravitational 
constant, M is the mass of the star, and r is the unit vec- 
tor in the radial dir ection. We omit the centrifugal force (cf. 
Kap yla et~aT1l201 lb|) . The rate of strain tensor S is given by 
Sij = | {ui-j +Uj.i) — ^5ij 1 V-u, where the semicolons denote 
co variant differentiation. 

2.1. Initial and boundary conditions 

In the initial state the atmosphere is adiabatic and the hy- 
drostatic temperature gradient is given by 
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where m = 1.5 is the polytropic index. We use Eq. © as 
the lower boundary condition for the temperature. The den- 
sity profile follows from hydrostatic equilibrium. The heat 
conduction profile is chosen so that radiative diffusion is re- 
sponsible for supplying the energy flux in the system, with K 
decrea sing more than two orders of magnitude from bottom 
to top dKapvla et al.ll2011ah . 

The radial and latitudinal boundaries are taken to be impen- 
etrable and stress free, according to 
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For the magnetic field we assume perfect conductors at the 
lower radial and latitudinal boundaries, and radial field at the 
outer radial boundary. In terms of A these translate to 
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On the latitudinal boundaries we assume that the thermody- 
namic quantities have zero first derivative, thus suppressing 
heat fluxes through the boundaries. 
On the upper boundary we apply a black body condition 
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where a is the Stefan-Boltzmann constant. We use a modi- 
fied value for a that takes into account that our Reynolds and 
Rayleigh numbers are much smaller than in reality, so K and 
therefore the flux are much larger than in the Sun. 



2.2. Dimensionless parameters 
We obtain non-dimensional quantities by choosing 

R = GM = po = cp = fiQ — 1 , 
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where po is the density at 0.7 R. The units of length, velocity, 
density, specific entropy, and magnetic field are then given by 



[x]=R, [u] = ^GM/R, [p} = 
[s] = c P , [B] = vWoGM/i?. 
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The simulation is governed by the Prandtl, Reynolds, Corio- 
lis, and (semi-) turbulent Rayleigh numbers, defined by 
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where Xtm is the turbulent thermal conductivity in the mid- 
dle of the convection zone (i.e. at r m = O.SbR), fcf = 
2n/Ar is an estimate of the wavenumber of the largest ed- 
dies, Ar = 0.3i? is the thickness of the layer, and ^ rms = 
\J (3/2)(u 2 , + Uq) is the rms velocity, where the angular 
brackets denote volume averaging. In our definition of u rms 
we omit the contribution from the azimuthal velocity, because 
its v alue is dominated by effects from the differential rota- 
tion ( Kap yla et al.ll201 lb|) . The magnetic field is expressed in 

equipartition field strengths, B eq (r) = {pvpu 2 ) 1 ^ where 
the subscripts indicate averaging over 6 and <fi, and az- 
imuthally averaged mean flows are subtracted. The Tay- 
lor number is Ta = (2Q J R 2 /u) 2 = Co 2 Re 2 (/c f iT) 4 , with 
kfR w 21. 

Due to the fact that the initial stratification is isentropic, we 
quote the value of Ra t from the thermally relaxed state of the 
run. The energy that is deposited into the domain at the base 
is controlled by the luminosity parameter 
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where Lq = 47rr 2 Fb is the constant luminosity, n = 0.7 R, 
and Fb = — (KdT/dr)\ r= o^R is the energy flux imposed at 
the bottom boundary. Furthermore, the stratification is deter- 
mined by the normalized pressure scale height at the surface 
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where Ty = T(r = R). Si milar parameter definitions were 
used bv lDobler eTaLl (|2006). We use f = 0.02 which results 
in a density contrast of 30. 

The simulations were performed with the PENCIL CODE 1 , 
which is a high-order finite difference method for solving the 
compressible equations of magnetohydrodynamics. 

3. RESULTS 

As a hydrodynamic progenitor run, we take a thermal ly re- 
laxed snapshot from Run B4 of Kapyla et~aTI d201 lah with 
Co = 7.6, Re = 36 (so Ta w 1.4 x 10 10 ), Ra t « 3 • 10 6 , 
Pr = 2.5, C = 3.8 • 10~ 5 and add a weak (10- 4 £ eq ) small- 
scale seed magnetic field. We use Pm = 1 and a resolution of 
128x256x128 mesh points. The magnetic field grows expo- 
nentially over roughly 1500 convective turnover times before 
reaching the saturated stage in which the total rms magnetic 
field is B rms = 0.72£ eq . 

The convection pattern near the surface (r = 0.98R) shows 
smaller scales at high latitudes and larger elongated struc- 
tures or 'banana cells' near the equator; see Fig. Q] Figure |2] 
shows the rotation profile in the saturated regime of the dy- 
namo. The equator rotates faster than the high latitudes and 
significant radial differential rotation is present only near the 
equator. In the lower part of the convection zone, 8ft / dr is 

1 http://code.google.eom/p/pencil-code/ 
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Fig. 1 . — Left: radial velocity u r normalised by the local sound speed c s 
near the surface of the star at r = 0.98R, duplicated fourfold in the azimuthal 
direction for visualization purposes. Right: radial velocity in the meridional 
plane for the same snapshot as in the left panel from <p = <f>o. 




Fig. 2. — Left: normalized time-averaged mean rotation profile Q/Qo = 
U / (flor sin 6) + 1 (contours) and meridional circulation U m = 
(U r ,Uo : 0) (arrows). Right: relative kinetic helicity h Te \. 

negative at low latitudes (just outside the inner tangent cylin- 
der) and positive at high latitudes. The meridional circulation 
shows several small cells outside the tangent cylinder on both 
hemispheres. The latitudinal differential rotation, measured 
by A n = (Qe=e - ^ eq )/^ eq , where ^ eq = Q(0 = tt/2), 
decreases from 0.08 in the kinematic regime to 0.07 in the 
saturated state. 

We define mean quantities as averages over longitude and 
denote them by an overbar. The relative kinetic helicity h re \ = 
u • ui/u rms uj rms , with uj = V x u, is negative (positive) 
in the northern (southern) hemisphere; see Fig. [2] No pro- 
nounced sign reversal with depth is seen. The maximum value 
of h re \ is around 0.3, which allows us to determine the dy- 
namo number describing the strength of the a effect as C a = 

a/Vtoh ~ h re \k^/ki « 2.7, where k^ = w rms /w rms is 
the approximate wavenumber of the energy-carrying eddies, 
fci = 7r/ Ar is the lowest radial wavenumber in the domain, 

while a « h re \u rms /3 and rjto = u rrns /3k^ are estimates 
for a effect and turbulent magnetic diffusivity. (We note that 
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Fig. 3. — Top panel: B^(r, t) in units of the local equipartition field 
strength at 25° latitude. Middle panel: near the surface of the star at 
r = 0.98R as a function of latitude 90° — 6. Bottom panel: blowup of the 
region -60° < 90° - < 60° and 2200 < tu rm sh < 3200. The white 
dotted line denotes the equator 90° — = 0. 

Xt / Vto varies between 1 .9 near the surface and 0. 15 within the 
convection zone.) The relevant dynamo number characteriz- 
ing the radial differential rotation is Cq = Aft/r] t okf w 55, 
where we have used AQ /Qq — 0.06 for the normalized radial 
shear (not to be confused with d efined above). The rati o 
Cn/Ca is well over 10. Following [Roberts & Stixl (119721) . 
this suggests that we are in what is known as the a£l regime 
where shear is strong enough to favour cyclic behavior. 

Time series of the averaged longitudinal component of the 
magnetic field are shown in Fig. Two activity belts are 
visible, one propagating poleward at high latitudes and an- 
other propagating equatorward between 10 and 30 degrees 
latitude. This mode becomes apparent in the nonlinear phase 
whereas in the kinematic stage the solution does not oscillate. 
Opposite transitions (oscillatory to quasi- steady) have been 
observed in Cartesian simulations ( Kapy lFet al 1 120 1 2h . The 
magnetic field is strongest at r/R w 0.85 and seems to prop- 
agate from there to top and bottom of the convection zone. 

Visualizations of the magnetic field near the surface (Fig.|4]) 
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t xi rms A, = 2140. 




Fig. 4. — Snapshots of the toroidal magnetic field at r = 0.98R from six different times separated by Atu rms kf « 105. 
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Fig. 5.— Top panel: (black line) and B r (red) at 90° - N = 25° 
latitude. The blue line shows 0.5£> r at #o- Bottom panel: B^ from #n and 
#S corresponding to latitudes ±25 degrees, respectively. 

show a persistent activity belt near the equator which is chang- 
ing polarity with a period of roughly Atu rms kf « 400. The 
same cyclic behaviour is seen throughout the depth of the 
convection zone. Relating the turnover time of our model, 
(w rms /cf) _1 , to that of the deep layers of the solar convec- 
tion zone, i.e. one month, leads to a magnetic cycle period of 
roughly 33 years. The cycle might well be shorter if the rel- 
evant depth is shallower. On the other hand, if we used k\ ^ 
instead of fcf , our cycle period would be 4-5 times longer. It is 
also noteworthy that has mixed parity about the equator, 



except around the time tu rms kf = 2500 when the field is of 
odd parity; see Fig. [5] 

On theoretical grounds, we would expect {B^/B^ to be 
of the order of IC^/Cy 1 / 2 « 4.5, but the actual ratio is 
only around unity; see Fig. [5] We cannot therefore be cer- 
tain that the dynamo is really in the aft regime, as discussed 
above. Moreover, B r shows a greater amplitude at high lat- 
itudes while B^ is stronger at lower latitudes. Furthermore, 
B r at high latitudes changes sign approximately when B^ in 
the low latitude activity belt changes sign. 

4. CONCLUSIONS 

We have reported solar-like magnetic cycles from a simula- 
tion of turbulent convection in spherical wedge geometry. The 
magnetic activity is concentrated in two belts, a high-latitude 
one propagating poleward, and a low-latitude one propagating 
equatorward. The strongest magnetic fields, however, occur 
in the high-latitude activity branch. Simulations with moder- 
ately slower and faster rotation show similar behavior. These 
results will be discussed in forthcoming publications. Relat- 
ing the convective turnover time in the simulation to that of 
the Sun we obtain a cycle period of 33 years which is some- 
what longer than t hat in the Sun and half that obtained by 
iGhizaru et al.l (|2010|) from quite a different model exhibiting 
similar solutions, but without equatorward migration. One of 
the differences to Kap yla et"aT] (120101) is the use of a black 
body boundary condition for the temperature. Furthermore, 
convective energy transport now dominates over radiative dif- 
fusion fcf. iKapyla et al.ll201 la|) . the density stratification in 
the convection zone is roughly an order of magnitude larger, 
and we have omitted a stably stratified overshoot layer be- 
neath. However, compared with the Sun, our contours of dif- 
ferential rotation are still too cylindrical and also the banana 
cell pattern of radial velocity might not be realistic. Both may 
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be a consequence of having a large Taylor number; even the 
turbulent Taylor number, 0/z/ t ) 2 Ta = 9Ta/Re 2 « 10 8 , 
is rather large. Here, v t « 7ft o has been used as an es- 
timate of the turbulent viscosity. The magnetic activity in 
our model is distributed throughout the convection zone, in 
contr ast to the widely accepted flux-t ransport dynamo mecha- 
nism (Dik pati & Charbonneaull999|) in which a one-cell anti- 
clockwise (north) meridional circulation is crucial. In our 
model, neither meridional flows, which are here too irregu- 
lar (see Figure 13, nor negative radial different ial rotation in 
the n ear- surface shear layer (as anticipated by iBrandenburgl 
120051) . which is here absent, seem to be able to explain the re- 
sulting equatorward migration. Clarifying this is an important 



goal for future work. 



The simulations were performed using the supercomputers 
hosted by CSC - IT Center for Science Ltd. in Espoo, Fin- 
land, who are administered by the Finnish Ministry of Educa- 
tion. Financial support from the Academy of Finland grants 
No. 136189, 140970 (PJK) and 218159, 141017 (MJM), as 
well as the Swedish Research Council grant 621-2007-4064, 
and the European Research Council under the AstroDyn Re- 
search Project 227952 are acknowledged. The authors thank 
NORDITA for hospitality during their visits. 



REFERENCES 



Brandenburg, A. 2005, ApJ, 625, 539 

Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 

2008, ApJ, 689, 1354 
— . 2010, ApJ, 711, 424 

Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 
2011, ApJ, 731, 69 

Browning, M. K., Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 648, 
L157 

Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 

Chan, K. L., & Sofia, S. 1986, ApJ, 307, 222 

Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 

Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 

Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 

Gilman, P. A. 1983, ApJS, 53, 243 

Glatzmaier, G. A. 1985, ApJ, 291, 300 

Kapyla, P. J. 2011, Astron. Nachr., 332, 43 

Kapyla, P. J., Korpi, M. J., Brandenburg, A., Mitra, D., & Tavakol, R. 2010, 

Astron. Nachr., 331, 73 
Kapyla, P. J., Korpi, M. J., & Tuominen, I. 2006, Astron. Nachr., 327, 884 



Kapyla, P. J., Mantere, M. J., & Brandenburg, A. 2011a, Astron. Nachr., 
332, 883 

— . 2012, submitted to GAFD, arXiv:l 11 1.6894 

Kapyla, P. J., Mantere, M. J., Guerrero, G., Brandenburg, A., & Chatterjee, 

P. 2011b, A&A, 531, A162 
Kitchatinov, L. L., & Olemskoy, S. V. 2011, Sol. Phys., 395 
Krause, F., & Radler, K.-H. 1980, Mean-field Magnetohydrodynamics and 

Dynamo Theory (Oxford: Pergamon Press) 
Miesch, M. S., Brown, B. P., Browning, M. K., Brun, A. S., & Toomre, J. 

2011, in IAU Symposium, Vol. 271, IAU Symposium, 261-269 
Miesch, M. S., & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317 
Ossendrijver, M. 2003, A&A Rev., 11, 287 
Parker, E. N. 1955, ApJ, 122, 293 

Racine, E., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, 

P. K. 2011, ApJ, 735,46 
Roberts, P. H., & Stix, M. 1972, A&A, 18, 453 



